{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 磁性质数值导数 (3)：RMP2 的 GIAO 核磁 (屏蔽) 共振常数 (NMR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 创建时间：2020-08-31；[修订前的文档](NMR_GIAO_NumDeriv_old.ipynb) 对原子轨道积分的磁化率梯度叙述上有错误。\n",
    ">\n",
    "> 修订时间：2022-02-09"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这篇文档中，我们会讨论使用 PySCF 以及其作为 libcint 的接口，计算 GIAO 的 RHF 数值核磁 (屏蔽) 共振常数 (Nuclear Magnetic Resonance constant, NMR) 的程序。该文档大量参考 PySCF 的代码 [nmr/rhf.py](https://github.com/pyscf/pyscf/blob/master/pyscf/prop/nmr/rhf.py)。\n",
    "\n",
    "该文档经过修订。之所以发现先前文档的错误，是因为读到最近的一篇数值 NMR 工作[^Glasbrenner-Ochsenfeld]。目前版本的文档应当能处理类似于 MP2、dRPA@HF、CCSD 等 Post-HF 方法。DFT 方法的 NMR 暂不在本文档的讨论范畴。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[^Glasbrenner-Ochsenfeld]: Glasbrenner, M.; Graf, D.; Ochsenfeld, C. Benchmarking the Accuracy of the Direct Random Phase Approximation and σ-Functionals for NMR Shieldings. *J. Chem. Theory Comput.* **2022**, *18* (1), 192–205. doi: [10.1021/acs.jctc.1c00866](https://doi.org/10.1021/acs.jctc.1c00866)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{warning}\n",
    "\n",
    "该文档已经修订完毕，但使用到了 [先前的一份文档](Mag_GIAO_NumDeriv.ipynb) 对 GIAO 的公式、记号与讨论。先前的文档还未修订完毕。\n",
    "\n",
    "之前犯的主要问题是，在解析导数求解时，将 Fock 贡献归入 Core Hamiltonian 是程序实现上非常方便的；但原理上，Fock 贡献要拆分为单电子与双电子部分。因此，在求数值导数时，自洽场部分的双电子积分 (ERIs) 也应是受了磁场影响的。\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "与之前的文档一样，我们的讨论中所使用到的分子体系 `mol` 会是非对称的氨分子，并且取用最小基组。其 RHF 计算放在实例 `mf`，而 NMR 计算实例会放在 `mf_nmr`。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 准备工作"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyscf import gto, scf, mp\n",
    "from pyscf.prop import nmr\n",
    "from pyscf.data import nist\n",
    "from scipy import constants\n",
    "import numpy as np\n",
    "\n",
    "np.set_printoptions(precision=5, linewidth=150, suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol = gto.Mole()\n",
    "mol.atom = \"\"\"\n",
    "N  0.  0.  0.\n",
    "H  0.  1.  0.2\n",
    "H  0.1 0.3 1.5\n",
    "H  0.9 0.4 -.2\n",
    "\"\"\"\n",
    "mol.basis = \"STO-3G\"\n",
    "mol.verbose = 0\n",
    "mol.build()\n",
    "coord_orig = np.zeros(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "nocc, nao, nmo, natm = mol.nelec[0], mol.nao, mol.nao, mol.natm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其自洽场能量为"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-55.253540514686556"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf = scf.RHF(mol).run()\n",
    "mf.e_tot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "RHF 的核磁屏蔽张量 (Shielding Constant) $\\sigma_{ts}^A$ 可以表示如下 (维度 $(A, t, s)$)：\n",
    "\n",
    "$$\n",
    "\\sigma_{ts}^A = \\frac{\\partial^2 E_\\mathrm{tot}}{\\partial \\mathscr{B}_t \\partial \\mu_{A_s}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[305.25641,  22.67382,  12.29714],\n",
       "        [-22.38193, 180.11491,  -2.84064],\n",
       "        [ 13.65836, -28.07439, 383.02872]],\n",
       "\n",
       "       [[ 25.68677,  -0.27932,  -0.31477],\n",
       "        [ -0.84315,  35.70697,  -4.26016],\n",
       "        [ -1.06334,   2.07184,  28.17535]],\n",
       "\n",
       "       [[ 24.06192,  -0.12066,   0.81476],\n",
       "        [  0.99821,  24.76595,  -0.86457],\n",
       "        [ -0.74954,   0.2574 ,  31.65924]],\n",
       "\n",
       "       [[ 39.48533,   4.28363,  -7.95068],\n",
       "        [  0.89957,  26.69818,  -0.9659 ],\n",
       "        [ -4.40617,  -0.48147,  28.90276]]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_nmr = nmr.RHF(mf)\n",
    "mf_nmr.kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上述结果是以 ppm 为单位的表达。不过我们在这篇文档打算考虑 RMP2 的核磁计算；这里的 RHF NMR 张量只是用于演示而已。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 微扰原子矩阵的表达"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 外磁场微扰的一阶 Core Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "关于该微扰矩阵，即是统合磁效应在 Core Hamiltonian 的作用、以及 GIAO 轨道对动能与核静电势能的一阶贡献考虑进来：\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "h_{\\mu \\nu}^{\\mathscr{B}_t}\n",
    "=\n",
    "\\frac{1}{2} \\langle \\mu | \\hat l_t | \\nu \\rangle_{\\mathrm{Gauge} \\rightarrow \\boldsymbol{R}_\\nu}\n",
    "+ \\langle U_\\mathrm{g}^t \\mu | \\hat t | \\nu \\rangle\n",
    "+ \\langle U_\\mathrm{g}^t \\mu | \\hat v_\\mathrm{nuc} | \\nu \\rangle\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "hcore_1_B = - 1j * (\n",
    "    + 0.5 * mol.intor('int1e_giao_irjxp', 3)\n",
    "    + mol.intor('int1e_ignuc', 3)\n",
    "    + mol.intor('int1e_igkin', 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 外磁场微扰的一阶重叠矩阵"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "S_{\\mu \\nu}^{\\mathscr{B}_t} = \\langle U_\\mathrm{g}^t \\mu | \\nu \\rangle\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们用下述代码生成 `ovlp_1_B` $S_{\\mu \\nu}^{\\mathscr{B}_t}$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "ovlp_1_B = - 1j * mol.intor(\"int1e_igovlp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 外磁场微扰的一阶 ERI 张量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(\\mu \\nu | \\kappa \\lambda)^{\\mathscr{B}_t} = (U_\\mathrm{g}^t \\mu \\nu | \\kappa \\lambda) + (\\mu \\nu | U_\\mathrm{g}^t \\kappa \\lambda)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们用下述代码生成 `eri_1_B` $(\\mu \\nu | \\kappa \\lambda)^{\\mathscr{B}_t}$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "eri_1_B = -1j * (\n",
    "    + np.einsum(\"tuvkl -> tuvkl\", mol.intor('int2e_ig1'))\n",
    "    + np.einsum(\"tkluv -> tuvkl\", mol.intor('int2e_ig1')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 核磁偶极的一阶 Core Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "核磁偶极 $\\mu_{A_t}$ 所产生的一阶算符贡献可以表达为\n",
    "\n",
    "$$\n",
    "\\hat h {}^{(1)} (\\boldsymbol{\\mu}_A)\n",
    "= - i \\alpha^2 \\boldsymbol{\\mu}_A \\cdot \\left( \\boldsymbol{\\nabla} \\frac{1}{\\boldsymbol{r - \\boldsymbol{R}_A}} \\times \\boldsymbol{\\nabla} \\right)\n",
    "= - i \\alpha^2 \\boldsymbol{\\mu}_A \\cdot \\left( \\frac{\\boldsymbol{r} - \\boldsymbol{R}_A}{|\\boldsymbol{r} - \\boldsymbol{R}_A|^3} \\times \\boldsymbol{\\nabla} \\right)\n",
    "$$\n",
    "\n",
    "其中，$\\boldsymbol{R}_A$ 表示原子 $A$ 的核坐标，$\\alpha$ 表示精细结构常数，$1/\\alpha \\simeq 137$。该常数可以从 PySCF 中获得，也可以从 SciPy 中获得。注意等式左边的 $\\boldsymbol{\\mu_A}$ 看作是外加的核磁偶极大小，而等号右边的 $\\mu$ 表示原子轨道，两者意义不同；等号左边角标 $\\boldsymbol{A}$ 表示原子核坐标向量，而之前两篇文档中的 $A$ 在很多文章或教材中表示 $\\frac{1}{2} \\boldsymbol{B} \\times \\boldsymbol{r}$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "137.0359990836958"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1 / constants.alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "但为了方便，在最后核算结果之前，我们会暂且将 $\\alpha$ 当作 1 来处理。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在 PySCF 中，实现上述过程的积分字符是 `int1e_ia01p`；但其使用需要告知 `gto.mole.intor` 函数以其 $1 / \\boldsymbol{r}$ 的规范原点位置具体处在哪个原子核中心。\n",
    "\n",
    "$$\n",
    "h_{\\mu \\nu}^{\\mu_{A_t}} = - i \\langle \\mu | \\boldsymbol{\\nabla} \\frac{1}{\\boldsymbol{r}} \\times \\boldsymbol{\\nabla} | \\nu \\rangle_{\\text{Gauge of } \\boldsymbol{r} \\rightarrow \\boldsymbol{R}_A}\n",
    "$$\n",
    "\n",
    "我们用下述代码生成 `hcore_1_m` $h_{\\mu \\nu}^{\\mu_{A_t}}$：(维度为 $(A, t, \\mu, \\nu)$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "hcore_1_m = np.zeros((natm, 3, nao, nao), dtype=np.complex128)\n",
    "for atom_idx in range(natm):\n",
    "    with mol.with_rinv_orig(mol.atom_coord(atom_idx)):\n",
    "        hcore_1_m[atom_idx] = - 1j * mol.intor(\"int1e_ia01p\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 磁场与核磁偶极的二阶 Core Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "磁场与核磁偶极之间的算符乘积会产生二阶算符贡献项：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat h {}^{(2)} (\\boldsymbol{\\mathscr{B}}, \\boldsymbol{\\mu}_A) | \\nu \\rangle\n",
    "&= \\frac{\\alpha^2}{2} \\boldsymbol{\\mathscr{B}}^\\mathrm{T} \\left( (\\boldsymbol{r} - \\boldsymbol{R}_\\nu) \\cdot \\boldsymbol{\\nabla} \\frac{1}{\\boldsymbol{r} - \\boldsymbol{R}_A} - (\\boldsymbol{r} - \\boldsymbol{R}_\\nu) \\boldsymbol{\\nabla} \\frac{1}{(\\boldsymbol{r} - \\boldsymbol{R}_A)^\\mathrm{T}} \\right) \\boldsymbol{\\mu}_A | \\nu \\rangle \\\\\n",
    "&\\quad + \\alpha^2 \\boldsymbol{\\mathscr{B}}^\\mathrm{T} \\boldsymbol{U}_\\mathrm{g} \\left( \\boldsymbol{\\nabla} \\frac{1}{\\boldsymbol{r} - \\boldsymbol{R}_A} \\times \\boldsymbol{\\hat p} \\right)^\\mathrm{T} \\boldsymbol{\\mu}_A | \\nu \\rangle\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "去除精细结构常数 $\\alpha$ 的贡献后，其矩阵的表达形式则是\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "h_{\\mu \\nu}^{A_s t} \\mathscr{B}_t \\mu_{A_s}\n",
    "&= \\langle \\mu | - \\frac{1}{2} \\frac{(t - t_\\nu) (s - s_A)}{|\\boldsymbol{r} - \\boldsymbol{R}_A|^3} | \\nu \\rangle\n",
    "- \\delta_{ts} \\langle \\mu | - \\frac{1}{2} \\sum_{w} \\frac{(w - w_\\nu) (w - w_A)}{|\\boldsymbol{r} - \\boldsymbol{R}_A|^3} | \\nu \\rangle \\\\\n",
    "&\\quad + \\langle U_\\mathrm{g}^t \\mu | \\left( \\boldsymbol{\\nabla} \\frac{1}{\\boldsymbol{r} - \\boldsymbol{R}_A} \\times \\boldsymbol{\\hat p} \\right)_s | \\nu \\rangle\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们用下述代码生成 `hcore_2` $h_{\\mu \\nu}^{A_s t}$：(维度为 $(A, t, s, \\mu, \\nu)$，注意维度 $t$ 对应外磁场，而 $A, s$ 对应核磁偶极)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "hcore_2 = np.zeros((natm, 3, 3, nao, nao))\n",
    "for atom_idx in range(natm):\n",
    "    with mol.with_rinv_origin(mol.atom_coord(atom_idx)):\n",
    "        hcore_2[atom_idx] += mol.intor(\"int1e_giao_a11part\").reshape((3, 3, nao, nao))\n",
    "        hcore_2[atom_idx] -= np.einsum(\"ts, uv -> tsuv\", np.eye(3), mol.intor(\"int1e_giao_a11part\").reshape((3, 3, nao, nao)).trace(axis1=0, axis2=1))\n",
    "        hcore_2[atom_idx] += mol.intor(\"int1e_a01gp\").reshape((3, 3, nao, nao))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数值导数求 RMP2 NMR 核磁屏蔽张量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "随后我们就可以通过数值梯度求核磁屏蔽张量 $\\sigma_{ts}^A$ (维度 $(A, t, s)$)：\n",
    "\n",
    "$$\n",
    "\\sigma_{ts}^A = \\frac{\\partial^2 E_\\mathrm{tot}}{\\partial \\mathscr{B}_t \\partial \\mu_{A_s}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在此之前，我们仍然需要构造一个通过更改 `get_hcore` Core Hamiltonian、`get_ovlp` 重叠矩阵、`_eri` 双电子积分的 PySCF 自洽场实例，以施加外场获得能量的函数 `eng_nmr_field`。其输入的参数 `dev_xyz_B` 是三维外加磁场大小 (对应维度 $t$)，`dev_xyz_m` 是三维外加核磁偶极大小 (对应维度 $s$)，`atom_idx` 是原子序号 (对应维度 $A$)。\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "h_{\\mu \\nu} &= h_{\\mu \\nu}^{(0)} + \\mathscr{B}_t h_{\\mu \\nu}^{\\mathscr{B}_t} + \\mu_{A_s} h_{\\mu \\nu}^{\\mu_{A_s}} + \\mathscr{B}_t \\mu_{A_s} h_{\\mu \\nu}^{\\mathscr{B}_t \\mu_{A_s}} + o(\\mathscr{B}_t \\mu_{A_s}) \\\\\n",
    "S_{\\mu \\nu} &= S_{\\mu \\nu}^{(0)} + \\mathscr{B}_t S_{\\mu \\nu}^{\\mathscr{B}_t} + o(\\mathscr{B}_t) \\\\\n",
    "(\\mu \\nu | \\kappa \\lambda) &= (\\mu \\nu | \\kappa \\lambda)^{(0)} + (\\mu \\nu | \\kappa \\lambda)^{\\mathscr{B}_t} + o(\\mathscr{B}_t)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在跑完自洽场后，立即简单地执行 MP2，就可以得到受外磁场与核磁偶极扰动的 MP2 能量了。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def eng_nmr_field(dev_xyz_B, dev_xyz_m, atom_idx):\n",
    "    mf = scf.RHF(mol)\n",
    "    def get_hcore(mol_=mol):\n",
    "        hcore_total  = np.asarray(scf.rhf.get_hcore(mol_), dtype=np.complex128)\n",
    "        hcore_total += np.einsum(\"tuv, t -> uv\", hcore_1_B, dev_xyz_B)\n",
    "        hcore_total += np.einsum(\"tuv, t -> uv\", hcore_1_m[atom_idx], dev_xyz_m)\n",
    "        hcore_total += np.einsum(\"tsuv, t, s -> uv\", hcore_2[atom_idx], dev_xyz_B, dev_xyz_m)\n",
    "        return hcore_total\n",
    "    def get_ovlp(mol_):\n",
    "        ovlp_total  = np.asarray(scf.rhf.get_ovlp(mol_), dtype=np.complex128)\n",
    "        ovlp_total += np.einsum(\"tuv, t -> uv\", ovlp_1_B, dev_xyz_B)\n",
    "        return ovlp_total\n",
    "    mf.get_hcore = get_hcore\n",
    "    mf.get_ovlp  = get_ovlp\n",
    "    mf._eri = mol.intor(\"int2e\") + np.einsum(\"tuvkl, t -> uvkl\", eri_1_B, dev_xyz_B)\n",
    "    mf.run()\n",
    "    mf_mp = mp.MP2(mf).run()\n",
    "    return mf_mp.e_tot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "随后使用与前文类似的数值梯度方式就能得到核磁屏蔽常数了；所使用的数值差分大小对于外磁场与核磁偶极均为 $3 \\times 10^{-4}$ a.u.。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "interval = 1e-4\n",
    "num_nmr = np.zeros((natm, 3, 3))\n",
    "for atom_idx in range(natm):\n",
    "    for t in range(3):\n",
    "        for s in range(3):\n",
    "            dev_xyzs_B, dev_xyzs_m = np.zeros((2, 3)), np.zeros((2, 3))\n",
    "            dev_xyzs_B[0, t] = dev_xyzs_m[0, s] = -interval\n",
    "            dev_xyzs_B[1, t] = dev_xyzs_m[1, s] =  interval\n",
    "            num_nmr[atom_idx, t, s] = (\n",
    "                + eng_nmr_field(dev_xyzs_B[0], dev_xyzs_m[0], atom_idx)\n",
    "                - eng_nmr_field(dev_xyzs_B[1], dev_xyzs_m[0], atom_idx)\n",
    "                - eng_nmr_field(dev_xyzs_B[0], dev_xyzs_m[1], atom_idx)\n",
    "                + eng_nmr_field(dev_xyzs_B[1], dev_xyzs_m[1], atom_idx)\n",
    "            ) / (4 * interval**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 5.63279,  0.41743,  0.26518],\n",
       "        [-0.11122,  3.81459, -0.13199],\n",
       "        [ 0.38918, -0.53118,  7.02831]],\n",
       "\n",
       "       [[ 0.47397, -0.00892,  0.00238],\n",
       "        [-0.02979,  0.68584, -0.06142],\n",
       "        [-0.0151 ,  0.0332 ,  0.52897]],\n",
       "\n",
       "       [[ 0.44192, -0.0033 ,  0.01888],\n",
       "        [ 0.0113 ,  0.45509, -0.01335],\n",
       "        [-0.01723,  0.00653,  0.59529]],\n",
       "\n",
       "       [[ 0.73128,  0.08291, -0.13904],\n",
       "        [ 0.03982,  0.50304, -0.0239 ],\n",
       "        [-0.06955, -0.01146,  0.54013]]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_nmr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "留意到我们之前一直都使用去除结构精细常数的结果，因此我们需要乘以 $\\alpha^2$。同时由于单位是 ppm，因此最终我们需要乘以 $10^6 \\alpha^2$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[299.95362,  22.2286 ,  14.12107],\n",
       "        [ -5.92259, 203.13221,  -7.02871],\n",
       "        [ 20.72447, -28.28604, 374.26729]],\n",
       "\n",
       "       [[ 25.2397 ,  -0.47482,   0.12685],\n",
       "        [ -1.58641,  36.52216,  -3.27057],\n",
       "        [ -0.80421,   1.76797,  28.16842]],\n",
       "\n",
       "       [[ 23.53299,  -0.17577,   1.00526],\n",
       "        [  0.60148,  24.23406,  -0.71087],\n",
       "        [ -0.91735,   0.3478 ,  31.69988]],\n",
       "\n",
       "       [[ 38.9417 ,   4.41492,  -7.40411],\n",
       "        [  2.12046,  26.78774,  -1.27294],\n",
       "        [ -3.70343,  -0.61013,  28.76271]]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_nmr_ppm = num_nmr * constants.alpha**2 * 10**6\n",
    "num_nmr_ppm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RMP2 NMR 核磁屏蔽常数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最终可以用于汇报的核磁屏蔽常数需要经过 $3 \\times 3$ 矩阵对角化给出。对于其中一个原子而言，其核磁屏蔽张量可以表示为\n",
    "\n",
    "$$\n",
    "\\boldsymbol{\\sigma} =\n",
    "\\begin{pmatrix}\n",
    "\\sigma_{xx} & \\sigma_{xy} & \\sigma_{xz} \\\\\n",
    "\\sigma_{yx} & \\sigma_{yy} & \\sigma_{yz} \\\\\n",
    "\\sigma_{zx} & \\sigma_{zy} & \\sigma_{zz}\n",
    "\\end{pmatrix}\n",
    "=\n",
    "X\n",
    "\\begin{pmatrix}\n",
    "\\sigma_{xx}' & 0 & 0 \\\\\n",
    "0 & \\sigma_{yy}' & 0 \\\\\n",
    "0 & 0 & \\sigma_{zz}'\n",
    "\\end{pmatrix}\n",
    "X^\\dagger\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中，$(\\sigma_{xx}', \\sigma_{yy}', \\sigma_{zz}')$ 是屏蔽张量的本征值，不随规范原点变化或坐标旋转而变化；而屏蔽张量 $\\boldsymbol{\\sigma}$ 本身是**随规范原点或坐标旋转变化而可以改变的**。因此，汇报 NMR 数据时，应当要汇报与本征值有关的结果。\n",
    "\n",
    "一般的 NMR 谱会打出同性核磁屏蔽常数 $\\sigma_\\text{iso} = \\frac{1}{3} (\\sigma_{xx}' + \\sigma_{yy}' + \\sigma_{zz}')$。有时异性屏蔽常数 $\\sigma_\\text{aniso} = \\sigma_{zz}' - \\frac{1}{2} (\\sigma_{xx}' + \\sigma_{yy}')$ 也会使用到[^Mason]。我们可以用下面的程序，对 MP2/STO-3G 的 NMR 屏蔽常数与{download}`Gaussian 结果 <NH3-nmr.out>` 作对照。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[^Mason]: Mason, J. Conventions for the Reporting of Nuclear Magnetic Shielding (or Shift) Tensors Suggested by Participants in the NATO ARW on NMR Shielding Constants at the University of Maryland, College Park, July 1992. *Solid State Nuclear Magnetic Resonance* **1993**, *2* (5), 285–288. doi: [10.1016/0926-2040(93)90010-K](https://doi.org/10.1016/0926-2040(93)90010-K)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      1  N    Isotropic =   292.5252   Anisotropy =   130.4969\r\n",
      "      2  H    Isotropic =    29.9780   Anisotropy =    10.0555\r\n",
      "      3  H    Isotropic =    26.4878   Anisotropy =     7.8257\r\n",
      "      4  H    Isotropic =    31.4998   Anisotropy =    15.9451\r\n"
     ]
    }
   ],
   "source": [
    "# Gaussian results\n",
    "! grep \"Isotropic\" NH3-nmr.out | tail -n 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def proc_nmr_tensor(tsr):\n",
    "    tsr = (tsr + tsr.T) / 2\n",
    "    eigs = np.linalg.eigvalsh(tsr)\n",
    "    eigs.sort()\n",
    "    nmr_iso = eigs.sum() / 3\n",
    "    nmr_ani = eigs[2] - (eigs[0] + eigs[1]) / 2\n",
    "    return nmr_iso, nmr_ani"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Isotropic = 292.4510  Anisotropy = 130.5998\n",
      "Isotropic =  29.9768  Anisotropy =  10.0491\n",
      "Isotropic =  26.4890  Anisotropy =   7.8233\n",
      "Isotropic =  31.4974  Anisotropy =  15.9435\n"
     ]
    }
   ],
   "source": [
    "# Numerical derivative results\n",
    "for tsr in num_nmr_ppm:\n",
    "    print(\"Isotropic = {:8.4f}  Anisotropy = {:8.4f}\".format(*proc_nmr_tensor(tsr)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
